/**
 * File: MinWeightedCoverageStrategy.java
 * Created by: mhaimel
 * Created on: 9 Jun 2009
 * CVS:  $Id: MinWeightedCoverageStrategy.java,v 1.2 2009/11/24 15:53:18 mhaimel Exp $
 */
package uk.ac.ebi.velvet.strategy;

import static uk.ac.ebi.curtain.utils.StatisticHelper.sum;
import static uk.ac.ebi.curtain.utils.StatisticHelper.wCov;

import java.util.ArrayList;
import java.util.Collections;
import java.util.List;

import uk.ac.ebi.velvet.processor.MinCutoffStrategy;

/**
 * @author mhaimel
 *
 */
public class MinWeightedCoverageStrategy implements MinCutoffStrategy {
	
	private static final double DEFAULT_STD = 0.1; // 10%
	private double std = 0.1;
	
	public MinWeightedCoverageStrategy() {
		this(DEFAULT_STD);
	}
	private MinWeightedCoverageStrategy(double std) {
		this.std = std;
	}

	/**
	 * {@inheritDoc}
	 */
	@Override
	public Integer calculateMinCutoff(MinimumCutoffInformation info) {
		int maxCov = info.getCoverageMaxLength();
		int sigma = getSigma(maxCov, info);
		int sigmaCov = maxCov-sigma;		
		
		List<Integer> covList = new ArrayList<Integer>(info.getCov2Score().keySet());
		Collections.sort(covList);
		
		Integer minCov = null;
		long minLen = 
			wCov(info.getCov2Score().get(maxCov),maxCov).longValue();
		
		for(int i = covList.size()-1; i >= 0; --i){
			int tmpCov = covList.get(i);
			if(tmpCov < maxCov && tmpCov < sigmaCov){
				long tmpLen = 
					wCov(info.getCov2Score().get(tmpCov),tmpCov).longValue();
				if(tmpLen < minLen){
					minCov = tmpCov;
					minLen  = tmpLen;
				}				
			}
		}
		return minCov;
	}
	
	private int getSigma(int maxCov, MinimumCutoffInformation info) {
		int iCnt = 0;
		double sig = 0.68;
		double sig2 = 0.95;
		long target = Math.round(info.getTotalLength()*sig);
		
		long temp = sum(info.getCov2Score().get(maxCov)).longValue();
		
		for(iCnt = 1; (maxCov - iCnt) > 1; ++iCnt){
			List<Integer> down = info.getCov2Score().get(maxCov-iCnt);
			List<Integer> up = info.getCov2Score().get(maxCov+iCnt);
			if(null != down){
				temp += sum(down).longValue();
			}
			if(null != up){
				temp += sum(up).longValue();
			}
			if(target <= temp){
				break;
			}
		}			
		int targetCov = maxCov - iCnt;
		if(targetCov > 1){
			// NOT targetCov
			return iCnt;
		} else {
			return (int)Math.round(Double.valueOf(maxCov*std));
		}
	}
}
